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The energetic optimization problem, e.g., searching for the optimal switching protocol of certain 
system parameters to minimize the input work, has been extensively studied by stochastic thermo- 
dynamics. In current work, we study this problem numerically with iterative dynamic programming. 
The model systems under investigation are toy actuators consisting of spring-linked beads with load- 
ing force imposed on both ending beads. For the simplest case, i.e., a one-spring actuator driven 
by tuning the stiffness of the spring, we compare the optimal control protocol of the stiffness for 
both the overdamped and the underdamped situations, and discuss how inertial effects alter the 
irreversibility of the driven process and thus modify the optimal protocol. Then, we study the sys- 
tems with multiple degrees of freedom by constructing oligomer actuators, in which the harmonic 
interaction between the two ending beads is tuned externally. With the same rated output work, 
actuators of different constructions demand different minimal input work, reflecting the influence of 
the internal degrees of freedom on the performance of the actuators. 

PACS numbers: 05.40.-a, 82.70.Dd, 87.15.H-, 05.70.Ln 
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I. INTRODUCTION 

The past two decades have seen the rapid growth of 
researches on nano-machines. Various types of nano- 
machines have been proposed and synthesized for prac- 
tical purposes [l|, This progress presents an urgent 
appeal and also a great challenge for physicists to un- 
derstand the working principle and energetics of nano- 
machines. For isothermal nano-machines including pro- 
tein machines and their derivatives or mimics, researchers 
are able to formulate a conceptual framework, the Brow- 
nian ratchet, to understand in general how a scalar en- 
ergy source (light, chemical reaction, thermal agitation, 
etc.) can facilitate a vectorial process. The ratchet 
mechanism has not only been used to explain many 
experimental observations on chemically-driven protein 
motors [1, 0, @, but also been successfully employed 
to develop some micro-manipulation techniques (e.g., a 
ratchet-like micro-device was designed for DNA segre- 
gation [6]). Within the ratchet framework, the ener- 
getics (or thermodynamics) has also been extensively 
investigated with multiple-state Langevin equations or 
Fokker-Planck equations (ej/., see [3[ for chemically- 
driven nano-machines, and p} for externally-controllable 
nano-machines) . 

In this context, some general energetic topics of nano- 
machines can be put forward. Particularly, in analogy 
to the optimization theory based on macroscopic finite- 
time thermodynamics [8, 9], the energetic optimization 
problems can be re-formulated for ratchet-type nano- 
machines. For instance, one may require the minimal 
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input work (abbreviated as MIW) in a finite-time pro- 
cess to achieve least heat agitation (lp| or maximum 
power [TH or highest efficiency [T^ |. To properly formu- 
late these questions, choosing proper model systems is 
very important. While the multiple-state ratchet model 
may not be a proper candidate (especially, not easy for 
analytical treatment), Sekimoto and Seifert et al. initi- 
ated the study with a much simpler ratchet model. In this 
model, the machine performs Brownian motion on an en- 
ergy landscape. And, a few parameters that characterize 
the landscape can be switched externally according to a 
deterministic protocol For such driven systems, 

stochastic thermodynamics has been firmly established, 
and the work and heat can be co-identified under the con- 
struction of first-law-like and second-law-like thermody- 
namic relations [HI, [l6|, [l7| • With the well-defined model 
systems and thermodynamics, a basic question is ad- 
dressed, i.e., what is the optimal protocol of the switching 
parameters to realize the minimal input work (abbrevi- 
ated as OPMIW). For the simple cases that the Brownian 
motion of an overdamped or underdamped bead is con- 
trolled by a tunable harmonic potential, the authors were 
able to present analytical solutions for the OPMIWs (see 
[ill [l3| for overdamped case, and [l4| for underdamped 
case). While their analysis focused on analytically solv- 
able potentials, Then and Engel also presented a Monte 
Carlo numerical method to discuss similar optimization 
problems of more complicated systems [loj . 

In this work, we follow the same logic and introduce a 
numerical method called iterative dynamic programming 
(abbreviated as IDP) 18] to investigate the OPMIW 
problems. For systems with polynomial potential func- 
tions up to the second order, the OPMIW problem is 
equivalent to the conventional optimal control problem, 
thus can be solved numerically with IDP. Inspired by 
the controllable DNA nano-actuator converting chemi- 
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cal energy into mechanical work jljj [2CJ , we choose toy 
actuators as our model systems. The simplest construc- 
tion is a two-bead actuator, in which the two beads are 
linked by a spring with externally tunable stiffness (a 
polymer can be regarded as such an actuator if its per- 
sistence length is tuned externally) . Although this model 
has been discussed by Seifert et al. [ll[ , the promising nu- 
merical method enables us to study the more complicated 
case when the actuator is working against a pulling force, 
and consequently to study the relation between the out- 
put work and the input work. Furthermore, since IDP 
can efficiently treat systems with multiple degrees of free- 
dom (l8j , actuators with different internal structures can 
be constructed and the related OPMIW problems can be 
discussed by this method. Besides, additional constraints 
can be easily specified (as will be clear in subsequent sec- 
tions) to study some interesting issues on energetics. 

In the following, we will first present the formulation of 
OPMIW problems for various toy models. In section [mi 
after testing the numerical precision of IDP for two ana- 
lytically solvable examples, we apply IDP to simple sys- 
tems with overdamped or underdamped dynamics, and 
also to systems with multiple degrees of freedom. Some 
general remarks will be presented in section IIVI The de- 
tails of IDP method can be found in the appendix. 



II. FORMULATION FOR OPMIW PROBLEM 

In this section, we will present the mathematical for- 
mulation of OPMIW problem for various model systems. 
We consider the simplest model of a one-spring actuator 
first, and then generalize the description to cases with 
multiple degrees of freedom. 

A. Formulation for one-spring actuator 

The one-spring actuator is basically the same model 
that Seifert and colleagues have considered before, in 
which a bead is trapped in a stiffening harmonic poten- 
tial [ll[. The actuator is illustrated in Fig. QJa), where 
two beads are connected by a spring with externally tun- 
able stiffness. To study the energetic problems, one bead 
is fixed at the origin while the other is pulled by a con- 
stant force. 

We only focus on one-dimensional problems here, thus 
the potential energy of the one-spring actuator can be 
written as 

V{x,X)= l -Xx 2 - fx, (1) 

where A is the stiffness of the spring, x is the position of 
the movable bead, and / is a constant pulling force (sup- 
pose / is pointing to the positive direction of z-axis). For 
the overdamped case, the Langevin equation describing 
the motion of the movable bead is 

7§ = -Az + / + v/2^T^(t), (2) 



where 7 is the frictional coefficient, ks is the Boltzmann 
constant, T is the temperature of the thermostat, and 
£(t) is Gaussian white noise satisfying 



(3) 



In the following, 7 and fcgT are set equal to unity for 
simplicity. If the unit of another quantity is specified, 
e.g., the force /, the dimensions of all the quantities in- 
volved in the model can be expressed by the dimensions 
of 7, fc_eT and /. The Fokker-Planck equation for the 
probability evolution of x reads 



dg(x,t) = d_ 
dt dx 



g(x,t) 



dV(x,X) , dg(x,t) 



dx 



dx 



(4) 



where g{x, t) is the probability distribution function of x 
at time t. 

According to stochastic thermodynamics (Tlj . under 
the framework of Langevin equation [Eq. p]) ]. the work 
increment dw along the stochastic trajectory of x is 
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aw = -^—Xdt = —Xx dt. 
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(5) 



Averaging with the Fokker-Planck equation [Eq. Q], 
when A is switched from an initial value Ai to a final 
value Xf within finite time, the expected input work Wi„ 
to the actuator can be expressed as 



Win = 



X(x 2 )dt, 



(6) 



where U and t / are the beginning and ending time point 
of the switching process respectively. The dot over A 
denotes the time derivative of A. The angular brackets 
denote an instant ensemble average using the probability 
distribution g(x,t). During the switching process, the 
average position (x) is changing, consequently the finite- 
time output work of the actuator is defined as 



W out = f(x)\t=t i -f(x)\ t= t r 



(7) 



There exist several differential constraints on the mo- 
ments of g{x,t) when optimizing Wt n . Multiplying both 
sides of Eq. (|4j) by x or x 2 and taking average, we obtain 
the following relations 



(x 2 ) = 2-2X{x 2 ) + 2f (x) 



(8) 



(9) 



In IDP calculation, we usually consider Wj„ as a function 
of time by defining 



W in {t)= J ^X(x 2 )dt, 
which satisfies the following differential equation 
Win = \u(x 2 ). 



(10) 



(11) 
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Here, the time derivative of A-protocol has been explicitly 
expressed as a temporal function u, i.e., 



X 



(12) 



Eq. ([nj and Eq. l[T2"|). together with Eq. © and Eq. J9|), 
compose the formulation of the OPMIW problem for cur- 
rent system. The task now turns into searching for the 
optimal trajectory of u that minimizes Wi n (tf). Once U, 
tf, Xi, Xf and the initial values of (x) and (x 2 ) are all 
specified [Wi„(0) is zero by definition], the problem can 
be solved with IDP method. 

For underdamped one-spring actuator, while Wi n is 
the same with Eq. ©, the set of differential constraints 
become 



(x 2 ) = 2(xp)/m, 



(13) 



(xp) = (p 2 )/m - (xp)/m - X(x 2 ) + f(x), (14) 



(p 2 ) = -2(p 2 )/m - 2X(xp) + 2f(p) + 2, (15) 



(x) = (p)/m, 



(p)=-(p)/m-X(x)+f. 



(16) 



(17) 



Where, p and m are the momentum and mass of the 
movable bead, respectively. 



Generalization to models with multiple degrees 
of freedom 



In this subsection, we consider a nano-system with co- 
ordinates vector x and time-dependent potential function 
V(x, X(t)). The parameter vector, A, of the potential is 
externally tuned according to certain protocol. To model 
the actuation process, constant forces are loaded on the 
system, thus the potential function can be written as 



V(x,X)=V (x, X)-f-x, 



(18) 



where, V$(x, A) represents the internal energy of the sys- 
tem (e.g., the elastic energy of the springs in the toy 
actuator) . / is not a normal three-dimensional force vec- 
tor, instead, it is a high-dimensional vector with each 
element representing the force on the corresponding de- 
gree of freedom. 

The dynamics of x can be described by overdamped 
or underdamped Langevin equations. For example, in 
overdamped case, the Langevin equations read 



1,2,.. 



,d, 



(19) 



where d is the dimension of x. For simplicity, we suppose 
that the random forces on different elements of x are 



uncorrelated, and the frictional coefficients are all the 
same. Thus, the noise terms (j = 1,2,..., d) satisfy 



(20) 



It is worth noting that models without such restrictions 
can also be solved with IDP. 

In the same spirit as section III Al when A is switched 
from its initial value A^ to the final value A/ in finite 
time, the input work to the system, Wi n , can be formally 
expressed as 



IT';,, - J X-(V x V(x,X))dt 



(21) 



where Va is the gradient with respect to A, and the brack- 
ets again denotes the average over the instant distribu- 
tion function g(x, t) of the system. Besides, during the 
switching, the change of (x) leads to an output work 



W oll 



f ■ (x)\ t =u ~ f ■ (x)\t=t f 



(22) 



If V(x, A) is a polynomial function up to the second or- 
der in x, the number of differential constraints on the 
moments of g(x, t) would be finite, the OPMIW problem 
can be rigorously expressed as an IDP-sovable optimiza- 
tion problem. 

To better illustrate the energetics of OPMIW, we only 
focus on the one-way switching processes of A throughout 
this article. One can integrate these and other processes 
into a complete working cycle to construct a novel nano- 
machine, as illustrated in Ref. [13] , the related optimiza- 
tion problem may need to be re-formulated. Besides, in 
present study, the systems are initially equilibrated un- 
der the potential function V(x,Xi). Consequently, the 
initial conditions of the differential constraints, i.e., the 
initial values of the moments of g(x,t), are selected as 
the corresponding equilibrium values. However, both the 
initial and final values of the moments can be chosen 
arbitrarily in computation as long as the optimal pro- 
tocol exists. Setting the final values of the moments of 
g(x, t) means that additional constraints are included in 
OPMIW problem. 



C. Formulation for oligomer actuators 

Under the general framework for systems with multi- 
ple degrees of freedom, we consider the one-dimensional 
oligomer actuators consisting of N+1(N — 1,2,...) beads 
in this subsection. Instead of the linear oligomers with- 
out branches and cross-links, we consider a sub-family of 
structured actuators with the potential function 

V(xi,x 2 , :.,x N ) = ^x\+^f =2 (x i -x i -x) 2 +-Xx 2 N -fx N . 

(23) 

Here, the beads are numbered from to N, Xi denotes 
the position of the ith bead. While the adjacent beads 
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are connected by a spring with elastic coefficient 1.0, xn 
is connected to Xq by an additional spring with tunable 
elastic coefficient A. To study the energetic problems, x 
is assumed to be fixed at the origin, a pulling force / 
is imposed on xn. The actuators are driven by switch- 
ing A from the initial value A^ to the final value A/ in 
finite-time. The current model system is ins pire d b y the 
recently designed DNA oligomer actuators [19|, |20|. In 
these DNA actuators, certain intra-molecular bonds can 
be modified by switching the light or chemical conditions, 
which results in considerable change of end-to-end dis- 
tance (i.e., a potential actuation process). Here, in the 
same spirit, a model actuator achieves an actuation when 
the 'intra-molecular bond' between bead and bead N 
is tuned externally. Therefore, the N = 1 model is ex- 
actly the above-mentioned one-spring actuator in which 
the only bond is tuned, and thus only the last two en- 
ergy terms in Eq. (|23[) need to be included in computa- 
tion. The N > 1 models describe structured actuators 
with extra degrees of freedom, e.g., the N = 2 actuator 
illustrated in Fig. QJb) . 

In current article, overdamped Langevin equations are 
chosen to model the dynamics of the actuators. For N — 
2 case, the optimal index, i.e., the input work, Wi n , is 

W in = f ' \x(x 2 2 )dt. (24) 

The differential constraints are 

(x\) = -4<z?) + 2{ Xl x 2 ) + 2, (25) 

(xix 2 ) = -(3 + X)(x 1 x 2 ) + (x\) + (xf) + f{xi), (26) 

(x\) = -2(1 + X)(xl) + 2(x x x 2 ) + 2f(x 2 ) + 2, (27) 

(a;'i) = -2(a;i) + (a; 2 ), (28) 

(x 2 )=-(l + X)(x 2 ) + (x 1 )+f. (29) 
For N — 3 case, Wi n is 

W in = j' \x{x\)dt. (30) 

The differential constraints are 

{x\)=-A{x\)+2(x 1 x 2 )+2, (31) 

(x\) = -A(x\) + 2( Xl x 2 ) + 2(x 2 x 3 ) + 2, (32) 

(xj) = -2(1 + X)(xl) + 2(x 2 x 3 ) + 2f(x 3 ) + 2, (33) 

(xix 2 ) = -A{X!X 2 ) + (xj) + (x\) + (jBix 3 ), (34) 



(£22:3) = -(3 + X)(x 2 x 3 ) + (xix 3 ) + (xf) + (x^) + f(x 2 ), 

(35) 

(xxx 3 ) = -(3 + X)(xix 3 ) + (x 2 x 3 ) + (xix 2 ) + f(xx), (36) 

(xi)=-2(xi) + (x 2 ), (37) 

(x 2 ) = -2(x 2 ) + ( Xl ) + (x 3 ), (38) 

(i 3 > = -(l + A)(a;3> + <a; 2 ) + /. (39) 

III. NUMERICAL RESULTS AND 
DISCUSSIONS 

A. Numerical accuracy of IDP 

To test the numerical accuracy of IDP, we first in- 
vestigate the two examples that have been analytically 
solved by Schmiedl and Seifert [TTj] . Both of the two 
systems are one-dimensional systems satisfying over- 
damped Langevin dynamics. The two model potentials 
are V\(x, A) = \(x— A) 2 (with tunable average position) 
and V 2 (x, A) = i Ax 2 (with tunable stiffness) respectively. 

The comparison between the numerical results and the 
analytical solutions is shown in Fig. [2l where the poten- 
tial functions, the initial and final values of A, the switch- 
ing time interval At and the MIW, Wi n , are all listed. We 
only list one MIW value in either figure, since the numer- 
ical MIW is the same with the analytical MIW [Eq. (9) 
and Eq. (19) in [111 ]] up to the shown numerical accuracy. 
In either figure, the solid line is the analytical protocol of 
A, and the hollow triangles denote the numerical results. 
Obviously, for the two examples studied here, the numer- 
ical results accurately follow the analytical curves. The 
important characteristics of the OPMIWs for the model 
systems, i.e., the jumps both at the beginning and at the 
end of the optimal protocols, are successfully reproduced. 

The triangles in Fig. [2] distribute non-uniformly along 
the time axis due to the employed method of variable 
step-length. In this scheme, the space between neighbor- 
ing time points are also optimized. Hence, more time 
points will be gathered at the sharply changing parts of 
the A-protocol, with less at the smoothly changing parts. 
As a result, higher accuracy can be obtained with limited 
discretization of [tj,i/]. Implementing the variable step- 
length method does not change the procedure of IDP (see 
the appendix). 

B. OPMIW for overdamped cases 

In this subsection, we consider the overdamped one- 
spring actuator introduced in section Hi Al Without the 
pulling force, this model turns into the analytical solvable 
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model with tunable stiffness (see section llll A [I . The force 
complicates the problem, making it very hard to get an 
analytical solution. However, the model system is quite 
suitable for IDP. For all the processes studied in this 
subsection, A is switched from 1.0 to 3.0, and the force is 
set as 1.0 if not further claimed. 

We first solve the OPMIW problems for processes with 
different switching time intervals and the same A^, A^- 
and /. Some of the OPMIWs are shown in Fig. [3][a). 
As can be seen, there are always jumps at the beginning 
and at the end of the optimal protocols, which is simi- 
lar to the cases in the already studied systems [13, E^l • 
Intuitively, with infinite or large enough time interval, 
all the jumps should approach zero no matter how large 
the force is (or even whatever the potential function is), 
because the switching process is actually quasi-static or 
reversible, 'reversible' here means that the input work 
can be completely converted into the free energy of the 
system, without any heat dissipation. This is consis- 
tent with the second law of stochastic thermodynamics 
Wi n > AF, where AF is the free energy difference be- 
tween the two equilibrium states that are determined by 
Xi and A / respectively. This inequality holds for any pro- 
tocols of control variables in a driven process including 
the OPMIW. To verify this, the MIWs for processes with 
different switching time are shown in Fig.[3jc). It can be 
seen that the MIW decreases monotonically with the in- 
crease of switching time, and asymptotically reaches AF 
when the switching time goes to infinity [the value of 
AF is indicated as the dashed line in Fig. [3jc)]. Con- 
sequently, the non-negative quantity Wi„ — AF can be 
taken as an indicator of the irreversibility of a driven pro- 
cess. In fact, it is exactly the total dissipation of the ac- 
tuator and the environment [la ] . In contrast to the input 
work, as shown in Fig. [3Jc), the output mechanical work 
increases monotonically with the switching time and also 
reaches an asymptotic value that is only determined by 
the initial and final equilibrium states. When the system 
finally reaches equilibrium, the change of average position 
(x) reaches its maximum with the largest output work. 
Thus, longer switching time results in a final state closer 
to equilibrium and a larger output work. The above re- 
sults suggest that IDP can provide physically sound and 
numerically correct results for the OPMIW problems dis- 
cussed here. 

While the infinitely slow process without any jumps 
leads to the lower bound (zero) of dissipation, the upper 
bound of dissipation, i.e., the highest irreversibility, oc- 
curs in the case of zero switching time (i.e., A is switched 
from Xi to Xf instantaneously). Between these two ex- 
tremes are the processes with finite switching time. Thus, 
it is reasonable to expect that the shorter switching time 
corresponds to the larger jumps in OPMIW [see Fig.[3fa)] 
due to the increased irreversibility. For switching pro- 
cesses with different time intervals, the relations between 
the magnitude of jumps (both of the initial and of the 
final jumps) and the irreversibility of the OPMIW are 
plotted in Fig. [3Jb). The trends of these curves indeed 



conform to the above guess. We also notice that when 
the switching time approaches zero, i.e., almost with the 
largest irreversibility, the OPMIW becomes stepwise and 
has very little advantage over the other control protocols: 
different A-protocols result in almost the same input work 
since the switching is too fast to change the position of 
the bead. 

Similarly, we also study the processes with the same 
switching time interval and various pulling forces. For 
current model system, larger force results in larger irre- 
versibility, which, as expected, leads to larger jumps in 
the optimized A-protocol [see Fig. (3Jd)]. From Fig. [3fb) 
and (d), together with other results of changing the 
switching range of A (data not shown) or the bead mass 
(the inertial effects will be discussed in the next subsec- 
tion), we expect that the magnitude of either the initial 
or the final jump in OPMIW may always monotonically 
correlate to the irreversibility, no matter what the poten- 
tial function is. Of course, this observation needs further 
validation in other systems. 



C. OPMIW for underdamped cases 

In this subsection, we study the underdamped one- 
spring actuator introduced in section ["II Al Similar to the 
overdamped cases, the OPMIWs for the underdamped 
actuators also have initial and final jumps. Besides, 
for the system with tunable average position in section 
IIII Al (5-function-like pulses have been analytically found 
in the optimized A-protocol. Although the current model 
system can not be solved analytically, the 5-function- 
like pulses in OPMIW have also been inferred when the 
pulling force is zero [lj| ■ For all the processes studied in 
this subsection, a constant force / = 1.0 is always loaded 
on the movable bead, A is switched from 1.0 to 3.0 with 
switching time 1.0. 

We first calculate the OPMIW for the one-spring ac- 
tuator with mass m = 0.02. Three optimal protocols 
with different allowable range of u = are plotted in 
Fig. 0Ja). In contrast to the overdamped cases, the op- 
timized A-protocols are no longer monotonic. There are 
always an upward peak at the beginning and a downward 
peak at the end of the OPMIW. The upward (downward) 
slope of either peak always matches the upper (lower) 
bound of u, indicating that the 5-function-like traces do 
exist in the OPMIW of current model system. The initial 
upward <5-peak enables an instant acceleration of the bead 
to ensure a nearly constant velocity during the process, 
and the final downward <5-peak enables an instant decel- 
eration of the bead to ensure minimal dissipation [Til ]. 
Nevertheless, to exactly reproduce a (^-function trace of 
A, the allowable range of u should be very large in com- 
putation, which is numerically hard (see the discussion 
in appendix) and is not the focus of this article. 

To demonstrate the feasibility of IDP-based OPMIW 
study for underdamped systems, we turn to a subset 
of protocols in which A grows monotonically with time. 
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This subset of protocols might be easier to implement in 
practice. We calculate the OPMIWs for actuators with 
different bead mass m, with two OPMIWs (to = 0.2 and 
to = 0.01) illustrated in Fig. OJa). The initial and fi- 
nal jump of A value are found again in the OPMIWs. 
As expected, when to decreases, thus j/m increases, the 
underdamped protocol will converge to the overdamped 
one. For finite to, however, plateaus appear immedi- 
ately after the initial jump and before the final jump 
in an optimized A-protocol. This phenomenon can be 
roughly explained as following. On one hand, after the 
first plateau stage, the bead has been accelerated to finite 
velocity. This velocity remains almost constant when A 
starts to increase again until reaching the other plateau 
of A-protocol. Therefore, from the perspective of minimal 
dissipation, the dissipation due to non-uniform velocity 
is minimized by including the first plateau stage. On the 
other hand, from the perspective of minimal input work, 
in the plateau stages with almost constant A, no work is 
input into the system. Besides, in the final plateau stage, 
the established velocity of the bead is damped. Thus, the 
extra dissipation after switching, which would be payed 
by the input work, is reduced. In a word, the protocol 
between the initial and final jumps can be well inter- 
preted either from the perspective of minimal dissipation 
or minimal input work. In the following, we will per- 
form more detailed analysis on the numerical OPMIWs 
to provide in-depth perspective of inertial effect, as well 
as to further verify the correctness of the computational 
results. 

As shown in Fig. [5fa), the OPMIW becomes stepwise 
when to increases. Taking the to = 0.2 actuator as ex- 
ample, except for the initial and final jump, a narrow 
region with sharply increasing A value exists approxi- 
mately in the middle of the OPMIW. Actually, this short 
piece of increasing A-protocol can be well approximated 
by a jump of A value without prominent increase of in- 
put work. In contrast, the position of this narrow region 
is more prone to changing the input work. To illustrate 
this point, we compare different stepwise protocols with 
three jumps and two plateau stages of A value. The ini- 
tial and final jumps of these protocols are set the same 
with the OPMIW of to = 0.2 actuator (thus the magni- 
tude of the intermediate jump is also determined), only 
the position of the intermediate jump is variable. For 
one-spring actuator with bead mass 0.2, the input works 
of these protocols are plotted in Fig. QJb). As can be 
seen, the stepwise protocol with the smallest input work 
has its intermediate jump located within the intermedi- 
ate increasing region of OPMIW, and its input work is 
only slightly larger than the MIW. This result strongly 
supports that the OPMIW obtained here is an optimum 
of the 'protocol sub-space' (i.e., the space of monoton- 
ical A-protocols). When to further increases, the input 
works of various A-protocols will gradually converge to 
the same value. This can be seen from Fig. [5jc), where 
the possible variation range of irreversibility is shown for 
actuators with different mass. The variation range is de- 



fined as the maximal irreversibility minus the minimal 
irreversibility (or equivalently, the maximal input work 
minus the minimal input work). The former is obtained 
by instantaneous switching and is obviously identical for 
actuators with different mass, and the latter is obtained 
with OPMIW. Fig. \5[c) shows that the larger bead mass 
corresponds to the larger minimal irreversibility and thus 
the smaller variation range, which means that all the A- 
protocols (including OPMIW) lead to nearly the same 
irreversibility and input work when to approaches infin- 
ity. 

The inertial effect can also be perceived from the rela- 
tion between input energy and system energy. For over- 
damped case, the potential energy of the system [i.e., 
the instant ensemble average of Eq. (Q])] increases mono- 
tonically with the cumulative input work as shown in 
Fig. [5]Jb). By contrast, as shown in Fig. [UJd), this rela- 
tion is no longer monotonic for the underdamped cases 
(only the to = 0.2 case is shown here). In the first plateau 
stage of OPMIW, while no work is input, part of the po- 
tential energy transforms into kinetic energy and heat, 
which is indicated by the first downward jump on the 
potential energy curve and the first upward jump on the 
kinetic energy curve in Fig. [5]Jd). During the ascend- 
ing stage in the middle of OPMIW, the kinetic energy 
remains almost constant, and the potential energy be- 
gins to increase with the cumulative input work. In the 
second plateau stage, both potential and kinetic energy 
are dissipated as heat, which is indicated by the second 
(downward) jump on the kinetic and potential energy 
curves. When to becomes larger, the second jump on ei- 
ther energy curve disappears first and then does the first 
jump (data not shown). When to approaches infinity, the 
bead is too heavy to be accelerated. Consequently, the 
whole system seems to be switched instantaneously from 
an equilibrium state to a non-equilibrium state. The in- 
put work is completely stored as potential energy and the 
kinetic energy remains unchanged. 



D. Oligomer actuator: OPMIW for more degrees 
of freedom 

In this subsection, we will investigate how the internal 
degrees of freedom can change the OPMIW and perfor- 
mance of an actuator. Concretely speaking, we try to 
study an optimal design problem, i.e., can we construct a 
more efficient actuator that output the specified amount 
of mechanical work with least input work? To explore 
this question, we study the model systems introduced in 
section III CI with N equals to I, 2 and 3 respectively. 
In a fast enough actuation process, an actuator might 
be driven out of equilibrium, thus its internal structure 
may begin to substantially affect the performance of the 
energetics, which leaves room for optimal design. 

Now, we can explicitly address the optimization prob- 
lem. With the same variation range of A € [A,-, A/], the 
same time interval t G [ij,t/] and identical pulling force, 



7 



as well as an extra requirement of equal amount of out- 
put work (i.e., rated output work), the MIW values for 
different structured actuators can be compared. The one 
with the lowest MIW should outweigh the others. For il- 
lustration, we only consider the overdamped N = 1,2,3 
models. In all the examples, A is switched from 1.0 to 
3.0, the switching time interval is [0.0,1.0], the pulling 
force is set as 1.0, and the rated output work is chosen 
as 0.3. 

The OPMIWs for the N = 1, 2, 3 models are plotted 
in Fig. [6ja). While there are both initial and final jumps 
in the OPMIWs of the N = 1 and N = 3 models, the 
final jump in the OPMIW of the N — 2 actuator disap- 
pears, and a plateau region appears instead. This is a 
consequence of constraining the output work. As shown 
in Fig.[6jb), among the three models considered here, the 
N = 2 actuator possesses a maximal finite-time output 
work which is the closest to the rated output work. Since 
this maximal output work is realized by instantaneously 
switching A from A^ to A/ at the beginning of the time 
interval (which is an infinitely fast process), the switch- 
ing process of the N = 2 actuator may also be relatively 
fast to realize the rated output work. Consequently, in 
contrast to protocols with final jump, A/ is reached be- 
fore t f in the OPMIW of N = 2 actuator. Actually, for 
any model system studied in this subsection, the OP- 
MIW will become one-step-like when the rated output 
work approaches the maximal finite-time output work. 

In spite of the fast switching OPMIW, the MIW of 
the N = 2 model turns out to be the smallest among 
the three models [sec Fig. GO^b)]. In other words, the 
N = 2 actuator performs the best when generating a 
finite-time output work of 0.3, which demonstrates the 
possibility for an optimal design of oligomer actuators. 
Interestingly, the N = 2 model also has the smallest irre- 
versibility among the three models [see the gap between 
MIW and free energy change in Fig.[Sfb)]. This observa- 
tion inspires us to propose an premature conjecture that 
irreversibility and performance might be related for ac- 
tuators with similar construction and different number 
of degrees of freedom. If such a relation does exist, the 
above discussion about optimal design could be gener- 
alized no matter the intra-molecular interaction is har- 
monic or not (for a real macromolecule, the intramolec- 
ular interaction is usually not harmonic), since it is the 
irreversibility that determines the energetic performance 
of an actuation device. 



at the end of OPMIW. For the underdamped case, the 
(5-function-like traces in OPMIW are also found. With 
vanishing mass, the underdamped models indeed con- 
verge to the overdamped one, and with increasing mass, 
the OPMIW of the underdamped models become step- 
wise. However, with quite large mass, the OPMIW only 
marginally outweigh the other protocols in terms of the 
input work, which is similar to other cases (such as van- 
ishing switching time) with large irreversibility. 

Furthermore, we use IDP to study the optimal con- 
trol and optimal design problems for the systems with 
multiple degrees of freedom. We consider a family of 
toy oligomer actuators in which the interaction between 
the two ending monomers can be tuned externally. The 
energetic performance of such actuators with different 
monomer number are compared, in the sense of out- 
putting the same amount of mechanical work with as- 
least-as-possible input work. Our results indicate that 
properly choosing the number of degrees of freedom can 
indeed bring better energetic performance. 

The above logic can be generalized to systems with 
more complex internal structures, e.g., a spring network 
with different harmonic interactions between the bead 
pairs. This construction is in the same spirit with the 
Gaussian Network Model (GNM) widely used to describe 
the thermal fluctuations of folded proteins [2l[ or the 
Elastic Network Model(ENM) used to describe the large- 
scale motion of motor proteins [22l ]. While those works 
focus on the structure-function and dynamics-function 
relations of protein machines, our studies might sug- 
gest a new perspective on energetics-structure relations. 
For instance, if the binding or release of substrates can 
be understood as tuning the local interaction between 
the residues of the motor proteins, energetic optimiza- 
tion problems similar to those studied here could also 
be raised. To carry out such investigation, however, 
all the thermodynamic relations and optimization prob- 
lems should be re-formulated for chemically-driven pro- 
tein machines. The stochastic thermodynamics of chemi- 
cal reaction networks [25| and the structural details of the 
macro-molecules should be combined, leading to a possi- 
bly formidable framework (one potential candidate is the 
generic ratchet model of chemically-driven machines Q) 
either for analytical or for numerical treatments. There- 
fore, at this preliminary stage, it may be beneficial to first 
construct some simple toy models to mimic the nano- 
machines. The toy actuators studied in this article can 
serve as the basic elements. 



IV. CONCLUSION AND PERSPECTIVE 

In current article, we study the OPMIW problems by 
introducing a numerical method, iterative dynamic pro- 
gramming (IDP). We first apply IDP to the systems with 
one degree of freedom and simple potential functions. 
Both the overdamped and the underdamped systems are 
investigated. For these two situations, jumps of the con- 
trol variable are always found both at the beginning and 



APPENDIX A: ITERATIVE DYNAMIC 
PROGRAMMING 

In the main text, we have presented the numerical re- 
sults of IDP for various OPMIW problems. Here we pro- 
vide a brief introduction to IDP. More details can be 
found in Ref. 

Generally speaking, IDP can solve the optimal control 
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problems of following style. 

f tf 

J = J L(s,u,t)dt, (Al) 



k=r(s,u,t). (A2) 

where L is a function of the state variables s, the con- 
trol variables u and time t. Since s are determined by 
u through the differential constraints Eq. (|A2[) . J is sim- 
ply an objective functional (also named as optimization 
index) of the protocol u(t). The related optimization 
problem is to find the optimal trajectory u(t) that min- 
imizes or maximizes J. Once the initial condition of the 
state variables and the initial time ti are specified, the 
problem is well formulated. In principle, it is also pos- 
sible to include more constraints on the final values of 
the state variables, the final time tf and the control vari- 
ables. Conventionally, this kind of problems is solved 
with Pontryagin principle (23|, though it is usually hard 
to obtain analytical solutions. The dynamic program- 
ming method (DP) [241] is a conceptually sound numeri- 
cal method for the optimal control problems. However, it 
also suffers from the difficulty of numerical interpolation 
and the restriction on dimensionality. In view of this, 
Luus developed IDP method to complement DP. With 
dramatic enhancement of speed, IDP is still promising 
for finding the global minimum, and achieves good nu- 
meri cal p recision as shown in current article and also in 
Ref. [18| . Here, we only describe the main procedures 
of IDP, readers can learn the underlying principles from 
Ref. [H. 

By introducing an additional state variable s* satisfy- 
ing 

s* = L{s,u,t), (A3) 
s*(U) = 0, (A4) 

s*(t f ) = J, (A5) 

the problem can be transformed into a set of differential 
equations shown in Eq. (|A6|) to Eq. (|A8|) . 

s* = L{s,u,t), (A6) 
s = r(s, u, t), (A7) 

J = «*(</). (A8) 

As a result, the optimization index J is determined solely 
by the state variables at time tf. 



1. The procedure of IDP 

The computational procedure of IDP is described in 
this subsection. 

Before calculation, the time interval, should be 

discretized into P stages with identical length I, i.e., 
PI = tf — ti, with constant control variables in each 
stage. Thus, the protocol of control variables can be 
represented by a P-dimensional array, each element of 
the array stores the values of control variables in the cor- 
responding time stage. In each time stage, an allowable 
range of control variables should be specified. The allow- 
able range does not change in one cycle of computation 
(the definition of a computation cycle will be given be- 
low). 

A reference P-dimensional array, U, of control vari- 
ables pre-exists to initiate each cycle. The elements of U 
are denoted as Ui (i = 1, 2, • • • , P). The allowable range 
of control variables is determined by U, together with the 
size of the range, R. For the case of one-dimensional con- 
trol, the allowable range in the ith time stage could be 
[Ui — R,Ui + R]. Of course, it is not necessary to set the 
reference control at the center of the allowable range. To 
start the first computation cycle, it is necessary to specify 
an initial P-dimcnsional array as the reference and also 
the initial size R of the allowable range. It should be 
noted that R could be either independent or dependent 
on time. We only discuss the former case here, since no 
substantial difference exists between the two cases. 

In calculation, a matrix of state variables should be 
first constructed as the mesh grid for dynamic program- 
ming. For simplicity, we only discuss the case of one- 
dimensional control in the following. 

First, The allowable interval of control variable at each 
time stage is uniformly divided into N s — 1 parts, with the 
N s nodes of this partition taken as the allowable values 
of control. All the largest allowable values in each time 
stage are grouped as the first P-dimensional array of con- 
trol (i.e., the first control protocol). The second largest 
ones are grouped as the second P-dimensional array, etc. 
Finally, N s P-dimensional arrays of control variable (N s 
control protocols) are obtained. Suppose the arrays are 
labeled with m, the elements of one array are labeled 
with n, then all these control values can be put into a 
N s x P matrix, with its elements denoted as u mn . 

Given the matrix of control variable, we can begin to 
construct the matrix of state variables. The idea is to use 
each P-dimensional array as the control protocol to cal- 
culate the corresponding state variables in different time 
stages. For the arbitrarily selected jth array, the calcula- 
tion starts from the initial state So, which is determined 
by the initial conditions. An updated state Sji can be 
obtained by taking Uj\ as the control variable and inte- 
grating Eq. |[A"6]) and Eq. (fAT]) from U to U + I. Sji is 
stored as an element into the matrix of state variables. 
After that, starting from Sji, another updated state Sj2 
can be obtained by taking Uji as the control variable and 
integrating from ti + 1 to ti + 21. With similar procedure, 
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the jth row of the matrix can be constructed by integrat- 
ing along the jth control protocol. After all the integra- 
tions are finished, we get a N s x P — 1-dimensional matrix 
of the state variables as shown in Fig. [7] Subsequently, 
the optimal control value for each element of the matrix 
can be estimated by the backward iteration method of 
dynamic programming. 

The backward iteration starts from the states 
Sip-i (i = 1,2, •• • ,N S ), i.e., the states at time £/ — I. 
T trial values of control variable should be provided in 
priori. Starting from certain Sip-x, taking respectively 
the T trial values as the control variable and integrating 
Eq. (|A"6|t and Eq. $K7} from tf - I to £/, we get T states, 
Sip-i U ' — 1j 2, • • • , T), at t = tf, which correspond to 
different J values. The trial control value leading to the 
minimal J (here we only consider the minimization prob- 
lems) is chosen as the optimal control variable for Sjp-i, 
and will be stored together with Sip—i. 

Subsequently for each Sip-2 (i = 1,2, ••• ,N S ), by 
one round of integration from tf — 21 to tf — I, T states, 
S\ P _ 2 (j = 1,2, • • • , T), at t = tf — I are generated. To 
proceed, we need to start from these states and integrate 
from tf — I to tf. Then, since the resulted T states are at 
t = tf, their corresponding J values can be compared to 
judge the optimal control value for S';p_2. Similarly with 
the preceding paragraph, the trial control value leading to 
the minimal J should be selected as the optimal control 
of SiP-2- The control values of S( P _ 2 (j = 1, 2, • • • , T) 
can be determined by the aid of the matrix of state vari- 
ables. The key is to find the nearest state to S 3 iP _ 2 among 

Sip-i (i = 1, 2, • • • , N s ). If the nearest state to S 3 P _ 2 is 

S m p-i, the control value of S 3 iP _ 2 should be set equal to 
the optimal control of S m p~\. 

The same procedure applies to other elements Sij (i = 
1, 2, • • • , N s .j = P - 3, P - 4, • • • , 1) of the matrix of 
state variables. One cycle is finished while the above 
procedure has been performed for So- After a cycle, 
we get the matrix with renewed optimal control val- 
ues u*j (i = 1, 2, • • • ,N s ;j = 1, 2, • • ■ , P — 1) for states 
Sij (i = 1,2, • • • , N s ; j = 1, 2, • • • , P — 1), as well as the 
optimal control value Uq for Sq. 

Before the next cycle, it is necessary to refresh the ref- 
erence array of control variable. The procedure is as fol- 
lowing. Starting from Sq, Eq. (|A6|) and Eq. (|A7j) are first 
integrated from ti to ti + I with the optimal control Uq, 
with a new state Sq obtained. Uq is stored as the first el- 
ement of a P-dimensional array. Then, within the states 
Sij (j — 1,2, ■•• , N s ), the one nearest to Sq is picked 
out. Suppose the nearest one is Si m , starting from Sq, 
and integrating Eq. (fA"6|) and Eq. (|a7|) from U+l to U+21 
with the optimal control value u\ m of Si m , we get a new 
state S^. Once again, u\ m is store in the P-dimensional 
array as the second element. A P-dimensional array of 
control variable can be generated by repeating the above 
procedure step by step. This protocol will be taken as 
the new reference array of control variable. Besides, the 
size of the allowable range, R, of control variable should 



be contracted before the new cycle. The new reference 
array and R together determine the allowable range of 
control in the next cycle. 

A pass is constituted by prescribed number of cycles. 
Several passes are included in each calculation. At the 
beginning of a new pass, R is restored partially to its 
initial value, and the optimal control protocol generated 
in last pass is used as the initial reference array of control 
variable. 



2. Variable step- length IDP 

The IDP algorithm with variable step-length is al- 
most the same with normal IDP method. Only slight 
re-formulation is needed. Here the time interval [£*,£/] is 
scaled to [0, 1]. As in normal IDP method, the [0, 1] inter- 
val is now uniformly divided into P segments. Further- 
more, each segment i corresponds to a real time length of 
Vi, and the real time length in each segment is introduced 
into IDP algorithm as an extra control variable. Due to 
the scaling of time, Eq. (|A6|) to Eq. (|A8|) for the ith. time 
segment should be expressed as below, 

%=L{s,u,t')v' i , (A9) 



— - r(s,u,t')— = r(WK, (A10) 

J = s*(l), (All) 

where, v[ satisfies 

<-£r«* < A12 > 

£' is the scaled time. The sum for all Vi should be equal 
to tf — ti. This constraint can be incorporated by sub- 
stituting the origin optimization index J by a new one I, 
which is the sum of J and a penalty function for violating 
the constraint. For example, / could be chosen as below 

1 = J + 6{Y li v i -t f +U) 2 . (A13) 

Where, 6 is a positive constant, Vi is the value of v at the 
ith time segment. The other kinds of constraints can be 
implemented in similar way. 

APPENDIX B: COMPUTATION DETAILS 

As pointed out in the main text, the time derivative of 
A-protocol is expressed as a temporal function u under 
current formulation. When searching for the OPMIW, 
u should be confined within a finite interval at all time. 
Consequently, larger range of u is necessary to describe 
the abrupt changing of A, e.g., a noncontinuous jump or 
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a (5-function trace of A value. However, if in a prelimi- 
nary calculation the optimized A-protocol is found to be 
monotonic, and some parts of the protocol are jump-like 
(for example, the u values in these parts are always equal 
to the allowable upper or lower bound), it is possible to 
accurately estimate the OPMIW in the following way. 
Instead of u = 4|, we can consider the reversed case, 
u' = 4j , in calculation because of the monotonicity of t 
versus A, and re-formulate all the functions of t to func- 
tions of A. Since there is usually no need to specify large 
range for u' , the calculation is simplified. This situa- 
tion happens to be true for all the overdamped examples 
studied in current article. Because the initial calculation 
without restriction almost ensures the monotonicity of 
A, the above-mentioned scheme should be promising for 
finding the global minimum. For demonstration, the re- 
formulated Eq. (dT]), Eq. dHJ), Eq. © and Eq. © of the 
overdamped one-spring actuator are shown below, 



dt 
dX 



d(x) 

~d)T 



fu — A (x) u , 



(Bl) 



(B2) 



(B3) 



In current article, we usually study each model (and 
parameter settings) with two IDP calculations. In the 
first calculation, the initial reference protocol is selected 
as a straight line connecting the two points, (tj,Aj) and 
(tf,Xf) (with constant u = Q or v! = After that, 
the resulted control protocol is selected as the initial ref- 
erence in the second calculation. According to our experi- 
ence, this procedure can ensure the successful application 
of IDP in all the examples here. 

In computation, we apply the Runge-Kutta method in 
Ref. [26j for integration. The error tolerance in integra- 
tion is set as l.e — 9 for the two examples in section ITlI A[ 
and l.e — 7 for the others. According to our experience, 
the MIW calculated with l.e — 9 tolerance are precise up 
to the seventh digit after the dot, the ones calculated with 
l.e — 7 tolerance are precise up to the fifth digit after the 
dot. Usually, the number of stage, P, is selected within 
the range of [15,20]. In the matrix of state variables, 
the number of state at each time point is selected within 
[10,20]. The number of trials along each dimension of 
control variable is about 10, leading to around 100 trials 
for two-dimensional control variables (including u = 4£ 
and the length of time steps). When new cycle begins, 
the contraction factor for the range of control variables is 
selected as 0.92. When new pass begins, the restoration 
factor for the range is selected as 0.88. Most of the cal- 
culations can be finished within 3 hours by single-thread 
computation with Intel(R) Q6700 CPU. 



—j- -L = 2u'-2\(x 2 )u' + 2f(x)u'. (B4) 
d\ 

It is worth noting that although we adopt the formulation 
with explicitly expressed time derivative of A, the other 
strategies that may be more natural for noncontinuous 
A-protocols, e.g., approximating the A-protocol with a 
stepwise function, can also be implemented with IDP. 
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(a) (b) 




FIG. 1: Illustration of the models in study, (a) The one-spring 
actuator with constant pulling force /. (b) The oligomer ac- 
tuator. The N — 2 case is portrayed. 
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FIG. 2: (color online). Comparison between numerical and 
analytical solutions [Eq.(9) and Eq.(19) in [ill ]] for over- 
damped systems with harmonic potentials, (a) Results for 
the system with tunable averaged position A. (b) Results for 
the system with tunable elastic constant (spring stiffness) A. 
The potential functions and related parameters are shown in 
the figure. 
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Time Irreversibility 



FIG. 3: (color online). Numerical results for the overdamped 
one-spring actuator, (a) The optimal protocols with switch- 
ing time 0.25 (squares), 1.00 (circles) and 10.0 (upward tri- 
angles). The switching time of the protocols have been re- 
scaled to [0.0, 1.0] for better visualization. The stiffness of 
the spring, A, is switched from 1.0 to 3.0. (b) Relation be- 
tween the initial (circles) or final (squares) jump in OPMIW 
and the irreversibility, Win — AF. The data points corre- 
spond respectively to processes with switching times of 25.0, 
10.0, 5.0, 3.0, 2.0, 1.5, 1.0, 0.75, 0.50, 0.25 and 0.10 (in the 
order of increasing irreversibility), (c) The minimal input 
works (squares) and the corresponding output works (circles) 
for processes with different switching time. The horizontal 
dashed line indicates the free energy difference between the 
two equilibrium states determined by Xi and A/. The hori- 
zontal dotted line indicates the maximal output work of the 
system, which can only be realized with infinite switching 
time, (d) Relation between the initial or final jump in OP- 
MIW and the irreversibility, Win — AF. The data points cor- 
respond respectively to processes with different pulling forces 
of 0, 0.25, 0.5, 1.0, 2.0, 5.0 and 10.0 (in the order of increasing 
irreversibility, the switching time is set as 1.0). 
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FIG. 4: (color online). The OPMIWs of underdamped one- 
spring actuator, (a) The non-monotonic OPMIWs for the 
m — 0.02 actuator with different bounds of u = \ ^\- The re- 
sults obtained with the bound of 5 (solid line), 10 (dashed 
line) and 15 (dash-dotted line) are shown. (b) The in- 
put works for protocols with three jumps and two plateaus 
(m — 0.2, see the main text for more details). The horizontal 
axis indicates the position of the intermediate jump. The left 
vertical axis shows the input work. The OPMIW (downward 
triangles) for this model is also shown for comparison. The 
right vertical axis shows the A value. The dotted horizontal 
line indicates the MIW value for this model. 




FIG. 5: (color online). Numerical results for the under- 
damped one-spring actuators, (a) The OPMIWs with differ- 
ent bead mass. The OPMIWs for m — 0.01 (hollow circles) 
and m = 0.2 (hollow triangles) are shown. For comparison, 
the OPMIW for the overdamped case is also shown (hollow 
squares). The curves have been shifted for better visualiza- 
tion. The relations between the instant energy terms of the 
system and the cumulative input work are shown in (b) (over- 
damped case) and (d) (m = 0.2), where the instant potential 
energy is denoted by hollow circles, the cumulative input work 
is denoted by dotted line, the kinetic energy is denoted by hol- 
low downward triangles. The horizontal axis is the cumulative 
input work. The arrows in (b) and (d) indicate the direction 
of time evolution, (c) shows the variation range of irreversibil- 
ity. The maximum and minimum irreversibility are plotted 
for different bead masses of 0.01, 0.5, 1.0, 5.0, 10.0 and 50.0. 
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FIG. 6: (color online). Numerical results for oligomer actu- 
ators, (a) The OPMIWs for the N = 1 (squares), N = 2 
(circles) and N = 3 (upward triangles) models. The switch- 
ing time is 1.0, and A is switched from 1.0 to 3.0. (b) The 
minimal input works (squares), the maximal output works 
in finite-time (upward triangles), the rated output work (cir- 
cles), and the free energy differences (downward triangles) for 
the N = 1, N = 2 and N = 3 actuators. 
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FIG. 7: The process for building the matrix of state variables. 



